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Abstract 

Hysteresis  in  smart  actuators  presents  a  challenge  in  control  of  these  actuators.  A  fundamental 
idea  to  cope  with  hysteresis  is  inverse  compensation.  In  this  paper  we  study  modeling,  identification 
and  inverse  control  of  hysteresis  in  smart  actuators  through  the  example  of  controlling  a  commer¬ 
cially  available  magnetostrictive  actuator.  The  (rate-independent)  Preisach  operator  has  been  used 
extensively  to  model  the  hysteresis  in  smart  actuators.  We  present  efficient  inversion  algorithms 
for  the  Preisach  operator  that  are  implementable  in  real-time.  The  magnetostrictive  hysteresis  is 
rate-dependent  at  high  frequencies.  For  this  we  propose  a  novel  dynamic  hysteresis  model  by  cou¬ 
pling  a  Preisach  operator  to  an  ordinary  differential  equation.  This  model  can  capture  the  dynamic 
and  hysteretic  behavior  of  the  magnetostrictive  actuator,  and  it  provides  insight  into  modeling 
of  rate-dependent  hysteresis  in  other  smart  materials.  The  effectiveness  of  the  identification  and 
inverse  control  schemes  is  demonstrated  through  extensive  experimental  results. 

Keywords:  Hysteresis;  Smart  actuators;  Preisach  operator;  Modeling;  Parameter  identification; 
Inverse  compensation 
‘Corresponding  author. 
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1  Introduction 


Smart  materials,  such  as  magnetostrictives,  piezoelectrics,  electroactive  polymers  (EAPs),  shape 
memory  alloys  (SMAs),  electrorheological  (ER)  fluids  and  nragnetorheological  (MR)  fluids,  all  dis¬ 
play  certain  coupling  phenomena  between  applied  electronragnetic/thernral  fields  and  their  mechani¬ 
cal/rheological  properties.  Actuators  and  sensors  made  of  these  materials  can  be  built  into  structures, 
often  called  smart  structures,  with  the  ability  to  sense  and  respond  to  environmental  changes  to  achieve 
desired  goals.  Smart  materials  and  smart  structures  have  been  receiving  tremendous  interest  in  the 
past  two  decades,  due  to  their  broad  applications  in  areas  of  aerospace,  manufacturing,  defense,  and 
civil  infrastructure  systems,  to  name  a  few.  Hysteresis  widely  existing  in  smart  materials,  however, 
makes  the  effective  use  of  smart  actuators  and  sensors  quite  challenging  [1],  In  this  paper  and  its 
sister  paper  [2],  we  study  modeling  and  control  of  hysteresis  in  smart  actuators,  through  the  example 
of  controlling  a  commercially  available  thin  magnetostrictive  actuator. 

Magnetostriction  is  the  phenomenon  of  strong  coupling  between  magnetic  properties  and  mechan¬ 
ical  properties  of  some  ferromagnetic  materials  (e.g.,  Terfenol-D):  strains  are  generated  in  response  to 
an  applied  magnetic  field,  while  conversely,  mechanical  stresses  in  the  materials  produce  measurable 
changes  in  magnetization.  This  phenomenon  can  be  used  for  actuation  and  sensing.  Magnetostrictive 
actuators  have  applications  to  micro-positioning,  robotics,  ultrasonics,  vibration  control,  etc.  Figure  1 
shows  a  sectional  view  of  a  Terfenol-D  actuator  manufactured  by  Etrerna  Products,  Inc.  By  varying 
the  current  in  the  coil,  we  vary  the  magnetic  field  in  the  Terfenol-D  rod  and  thus  control  the  motion 
of  the  rod  head.  To  give  some  idea  about  the  performance  of  magnetostrictive  actuators,  we  list  some 
specifications  for  the  AA-050H  series  actuators  manufactured  by  Etrerna  [3]:  Terfenol-D  rod  length  50 
mm,  maximum  dynamic  force  ±  220  N,  displacement  range  50  /in r.  Figure  2  displays  the  hysteresis 
observed  in  the  magnetostrictive  actuator. 

A  fundamental  idea  in  coping  with  hysteresis  is  to  formulate  the  mathematical  model  of  hysteresis 
and  use  inverse  compensation  to  cancel  out  the  hysteretic  effect.  This  idea  can  be  found  in  [5,  6,  7,  8, 
9,  10].  There  have  been  a  few  monographs  devoted  to  modeling  of  hysteresis  and  study  of  dynamical 
systems  with  hysteresis  [11,  12,  13,  14,  15]. 

Hysteresis  models  can  be  roughly  classified  into  physics-based  models  and  phenomenological  mod¬ 
els.  An  example  of  a  physics-based  model  is  the  Jiles- Atherton  model  of  ferromagnetic  hysteresis  [16], 
where  hysteresis  is  considered  to  arise  from  pinning  of  domain  walls  on  defect  sites.  The  most  popular 
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Stainless  Steel  Push  Rod 


Figure  1:  Sectional  view  of  a  Terfenol-D  actuator  [4] (Original  source:  Etrerna  Products,  Inc.). 


Figure  2:  Hysteresis  in  a  magnetostrictive  actuator. 

phenomenological  hysteresis  model  used  in  control  of  smart  actuators  has  been  the  Preisach  model 
[17,  5,  18,  19,  20,  9,  10].  A  similar  type  of  operator,  called  Krasnosel’skii-Pokrovskii  (KP)  operator  has 
also  been  used  [21,  8].  Although  in  general  the  Preisach  model  does  not  provide  physical  insight  into 
the  problem,  it  provides  a  means  of  developing  phenomenological  models  that  are  capable  of  producing 
behaviors  similar  to  those  of  physical  systems  (see  Mayergoyz  [12]  for  an  excellent  exposition). 

In  this  paper  we  study  modeling,  identification  and  inverse  control  methods  for  smart  actuators. 
We  have  chosen  a  magnetostrictive  actuator  as  an  example  for  two  purposes: 

The  magnetostrictive  hysteresis  is  rate-independent  when  the  input  frequency  is  low  (typically 
below  5  Hz):  roughly  speaking,  the  shape  of  the  hysteresis  loop  is  independent  of  the  rate  of 
input  variation  (see  Section  2  for  a  precise  statement).  The  rate-independent  hysteresis  can  be 
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modeled  by  a  Preisach  operator.  Note  that  a  variety  of  other  smart  actuators  have  been  modeled 
by  essentially  a  Preisach  operator  alone,  see,  e.g.,  [5,  20].  Therefore  the  identification  and  control 
methods  presented  here  apply  directly  to  a  wide  class  of  smart  actuators. 

-  The  magnetostrictive  hysteresis  is  rate- dependent  1  when  the  input  frequency  gets  high,  and  it 
can  no  longer  be  modeled  by  a  Preisach  operator  alone.  In  this  case  we  propose  a  novel  dynamic 
hysteresis  model,  consisting  of  a  Preisach  operator  coupled  to  an  ordinary  differential  equation 
in  an  unusual  way.  This  model  captures  important  dynamic  effects  in  the  frequency  region  of 
practical  interest.  We  expect  our  studies  on  rate-dependent  magnetostrictive  hysteresis  to  shed 
light  on  modeling  and  control  of  rate-dependent  hysteresis  in  other  smart  materials. 

The  remainder  of  the  paper  is  organized  as  follows.  We  provide  an  introduction  to  the  Preisach 
operator  in  Section  2.  Section  3  is  devoted  to  the  rate-independent  case.  We  first  discuss  an  iden¬ 
tification  scheme  in  Subsection  3.1,  and  then  present  inversion  algorithms  for  a  discretized  Preisach 
operator  and  for  a  Preisach  operator  with  nonsingular  measure  in  Subsections  3.2  and  3.3,  respectively. 
Modeling  and  control  of  rate-dependent  hysteresis  is  addressed  in  Section  4.  A  dynamic  hysteresis 
model  is  proposed  in  Subsection  4.1,  a  parameter  identification  method  is  presented  in  Subsection  4.2, 
and  an  inversion  algorithm  is  given  in  Subsection  4.3.  Finally  concluding  remarks  are  provided  in 
Section  5. 

A  preliminary  version  of  some  results  in  this  paper  was  presented  in  [22] . 


2  The  Preisach  Operator 


For  a  pair  of  thresholds  ((3,  a)  with  (3  <  a,  consider  a  simple  hysteretic  element  7 /3,a[',  •],  as  illustrated 
in  Figure  3.  For  u  €  C([0,T])  and  an  initial  configuration  (  £  {—1, 1},  the  function 


v  =  7/3, aKC]  :  [0,T]  ->  {-1,1} 


is  defined  as  follows  [13]: 


u(0)  =  <( 


-1 


if  u( 0)  <  (3 
if  (3  <  w(0)  <  a  j 


y  1  if  u(0)  >  a 

1In  some  literature,  e.g.,  [13,  14],  the  word  hysteresis  is  referred  to  rate-independent  memory  effects  only.  We  use 


“hysteresis”  in  a  more  general  sense  in  this  paper  and  [2]. 


4 


and  for  t  G  (0,  T],  setting  Xt  =  {r  E  (0 ,t]  :  u(t)  =  /3  or  a}, 


,  .  A  , 
«(*)  =  < 


«(0) 

-1 

1 


if  Xt  =  0 

if  X*  7^  0  and  it(maxXt)  =  f3 
if  Xt  /  0  and  ii(maxlt)  =  a 


v 


+1 

p 

a 

"  -1 

Figure  3:  The  elementary  Preisach  hysteron. 


This  operator  is  sometimes  referred  to  as  an  elementary  Preisach  hysteron  (we  will  call  it  a  hysteron 
in  this  paper),  since  it  is  a  building  block  for  the  Preisach  operator. 

The  Preisach  operator  is  a  weighted  superposition  of  all  possible  hysterons.  Define 

V0  =  {(/?,«)  E  M2  :  (3  <  a}. 

Vo  is  called  the  Preisach  plane ,  and  each  (f3,a)  G  Vo  is  identified  with  the  hysteron  7^.  For  u  G 
C([0,  T])  and  a  Borel  measurable  initial  configuration  ()o  of  all  hysterons,  (,0  '-Vo  {—1, 1},  the  output 
of  the  Preisach  operator  T  is  defined  as  [13]: 

z(t)  =T[u,(o\{t)  =  [  Ji3,a[u,(o{P,ai)}{t)dv(/3,a),  (1) 

JVo 

where  u  is  a  finite  (signed)  Borel  measure  on  Vo,  called  the  Preisach  measure. 

In  this  paper  we  call  the  Preisach  measure  v  nonsingular  if  \v\  is  absolutely  continuous  with  respect 
to  the  two-dimensional  Lebesgue  measure.  By  the  Radon-Nikodym  theorem  [23],  if  u  is  nonsingular, 
there  exists  a  Borel  measurable  function  /a,  such  that 

rKCo ]{t)=[  [  in{P,a)*/pta[u,(o{P,a)}{t)d(3da.  (2) 

J  Jv0 

The  weighting  function  p  is  often  referred  to  as  the  Preisach  function  [12]  or  the  density  function  [14]. 

To  simplify  the  discussion,  throughout  the  paper  we  assume  that  v(E)  =  0  for  any  Borel  measurable 
set 

E  C  {{j3,  a)  G  Vo\/3  <  /?o  or  a  >  oq} 
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for  some  /?o,  ao-  Then  it  suffices  to  consider  a  finite  triangular  area  V  =  {(/?,  a)  G  Vo\(3  >  An  a  <  07}, 
as  shown  in  Figure  4(a). 

The  memory  effect  of  the  Preisach  operator  can  be  captured  by  the  memory  curves  in  V.  At  time 
t,  V  can  be  divided  into  two  regions: 

V+(t)  =  {(/?,  ck)  G  V\  output  of  7 p>a  at  t  is  +  1}, 

V-(t)  =  {(/?,  ck)  G  V\  output  of  at  t  is  —  1}. 

Now  assume  that  at  some  initial  time  to,  the  input  u(to)  =  uq  <  /?o -  Then  the  output  of  every 
hysteron  is  —1.  Therefore  V-{to)  =  V,  V+(to)  =  0  and  this  corresponds  to  the  “negative  saturation” 
(Figure  4(b)).  Next  we  assume  that  the  input  is  monotonically  increased  to  some  maximum  value  at 
t[  with  u(ti)  =  u\.  The  output  of  7 pi0l  is  switched  to  +1  as  the  input  u(t)  increases  past  a.  Thus  at 
time  1 1,  the  boundary  between  V-{t\)  and  V+(ti)  is  the  horizontal  line  a  =  u\  (Figure  4(c)).  Next  we 
assume  that  the  input  starts  to  decrease  monotonically  until  it  stops  at  t,2  with  u(t 2)  =  U2 ■  It’s  easy 
to  see  that  the  output  of  7,3 )Cl  becomes  —1  as  u(t)  sweeps  past  (3,  and  correspondingly,  a  vertical  line 
segment  (3  =  U2  is  generated  as  part  of  the  boundary  (Figure  4(d)).  Further  input  reversals  generate 
additional  horizontal  or  vertical  boundary  segments. 


Figure  4:  Memory  curves  in  the  Preisach  plane. 

From  the  above  illustration,  each  of  V _  and  V+  is  a  connected  set  [12],  and  the  output  of  F  is 
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determined  by  the  boundary  between  and  V+  if  the  Preisach  measure  is  nonsingular.  The  boundary 
is  called  the  memory  curve.  The  memory  curve  has  a  staircase  structure  and  its  intersection  with  the 
line  a  =  (3  gives  the  current  input  value.  The  memory  curve  ifo  at  t  =  0  is  called  the  initial  memory 
curve  and  it  represents  the  initial  condition  of  the  Preisach  operator.  We  denote  by  \k  the  set  of  all 
memory  curves.  For  a  precise  characterization  of  'h,  we  refer  to  [24]. 

If  the  Preisach  measure  is  nonsingular,  we  can  identify  a  configuration  of  hysterons  Qi,  with  a 
memory  curve  if  in  the  following  way:  <0/,(/3,  a)  =  1  (— l,resp.)  if  (/?,  a)  is  below  (above,  resp.)  the 
graph  of  i ft.  Note  that  it  does  not  matter  whether  (0,  takes  1  or  —1  on  the  graph  of  if.  In  the  sequel 
we  will  put  the  initial  memory  curve  ifo  as  the  second  argument  of  T,  where  T[-,  ipo]  =  r[-,  <0O]. 

Theorem  2.1  summarizes  some  basic  properties  of  the  Preisach  operator,  see,  e.g.,  [13]. 

Theorem  2.1  [13]:  Let  v  be  a  Preisach  measure.  Let  u,u\,U2  £  C([0,T])  and  ifo  £  'k.  Then  the 
following  hold: 

1.  (Rate  Independence)  If  cf  :  [0,  T]  — >  [0,  T]  is  an  increasing  continuous  function  satisfying 

0(0)  =  0  and  0(T)  =  T, 

then  r[«  o  0,  ifo\(t)  =  r[ii,'0o](0(^)))  Vt  £  [0,T],  where  “o”  denotes  composition  of  functions. 

2.  (Strong  Continuity)  If  v  is  nonsingular,  then  r[-,0o]  :  C([0,T])  — >  C([0,T])  is  strongly 
continuous  (in  the  sup  norm). 

3.  (Piecewise  Monotonicity)  Let  v  >  0.  If  u  is  either  nondecreasing  or  noninreasing  on  some 
interval  in  [0,T],  then  so  is  r[u,0o]- 

4-  (Order  Preservation)  Let  v  >  0.  If  u\  <  U2  on  [0,  T],  then  r[wi,0o]  <  P['“2,  V;o]  on  [0,  T]. 


3  Identification  and  Inversion  of  the  Preisach  Operator 

3.1  Identification  of  the  Preisach  operator 

The  only  parameter  of  a  Preisach  operator  is  the  Preisach  measure.  A  classical  method  for  measure 
identification  is  to  use  the  so  called  first  order  reversal  curves  ,  detailed  in  Mayergoyz  [12].  Since 
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Figure  5:  (a)  Discretization  of  the  (restricted)  Preisach  plane  ( L  =  3);  (b)  A  memory  curve  (repre¬ 
sented  by  the  bolded  lines). 

the  method  involves  twice  differentiation,  a  smooth  approximating  surface  is  fit  to  the  data  points  in 
practice  [5,  18,  20].  As  pointed  out  in  [20],  deriving  the  density  by  differentiating  a  fitted  surface  is 
inherently  imprecise,  since  different  types  of  approximating  functions  lead  to  quite  different  density 
distributions.  Hoffmann  and  his  colleagues  proposed  a  scheme  to  identify  the  Preisach  measure  masses 
for  a  discretized  Preisach  operator  directly  [25] ,  based  on  which  the  density  function  is  approximated 
in  terms  of  a  set  of  basis  functions  [26] .  Another  identification  method  is  to  drive  the  Preisach  operator 
with  a  reasonably  “rich”  input  and  then  determine  the  measure  by  a  least  squares  algorithm  [21,  8,  9]. 
Here  we  briefly  review  the  identification  scheme  in  [9]. 

Smart  actuators,  due  to  the  capacity  of  the  windings  or  other  practical  reasons,  have  to  be  operated 
with  their  inputs  within  specific  ranges.  As  a  consequence,  one  can  only  visit  and  identify  the  measure 
distribution  in  a  restricted  area  A  in  the  Preisach  plane.  We  discretize  the  input  range  [umin,umax] 
into  L  +  1  levels  uniformly,  and  call  this  discretization  of  level  L  (see  Figure  5(a)  for  an  example  of 
L  =  3) .  The  Preisach  measure  within  each  cell  is  assumed  to  concentrate  as  a  discrete  mass  at  the  cell 
center  (see  the  dark  dots  in  Figure  5(a)).  This  corresponds  to  a  discretized  Preisach  operator,  which 
is  weighted  sum  of  a  finite  number  of  hysterons,  as  illustrated  in  Figure  6. 

To  identify  the  measure  masses,  one  applies  a  “rich”  input  function  u{t)  and  measures  the  output 
trajectory  z{t).  An  input  u{t)  is  rich  if  it  is  able  to  single  out  the  contribution  of  each  hysteron  in 
the  discretized  Preisach  operator,  and  one  candidate  for  such  u(t)  is  the  concatenation  of  the  first 
order  reversal  inputs.  Recall  that  a  first  order  reversal  input  is  obtained  by  increasing  the  input 
from  Umin  to  some  a  monotonically  and  then  decreasing  it  to  some  f3.  Signals  u(t),  z(t )  are  sampled 
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Figure  6:  The  discretized  Preisach  operator. 

into  sequences  {tt[n]}^=1,  {z[n\}^=l.  The  input  sequence  {w[n]}  (after  discretization)  is  fed  into  the 
discretized  Preisach  operator,  and  the  state  of  each  hysteron,  {7fc[n]},  k  =  1,  •  •  •  ,  K,  is  then  computed, 
where  K  is  the  number  of  hysterons.  The  output  of  the  Preisach  operator  at  time  instant  n  can  be 
expressed  as: 

K 

z[n\  =  r[u[-],^o]N  =  vo  +  y^*VXfcM;,  (3) 

k= 1 

where  vq  is  the  constant  contribution  of  the  Preisach  measure  from  outside  of  A,  and  v\~  with  k  >  0 
is  the  measure  mass  corresponding  to  %. 

Remark  3.1  We  put  a  sequence  instead  of  a  continuous  time  function  as  the  first  argument  ofT  in 
(3).  To  avoid  ambiguity,  it  is  tacitly  understood  that  the  input  is  changed  monotonically  from  u[n ] 
to  u[n  +  1].  Throughout  the  paper  we  may  use  a  sequence  or  a  continuous  time  function  as  the  first 
argument  ofT  depending  on  the  context. 


We  use  the  least  squares  method  to  estimate  the  parameters,  i.e.,  the  parameters  are  determined 
in  such  a  way  that 


N 

Ei 

n=l 


z\n\  —  z\n\ 


(4) 


is  minimized.  Since  we  require  >  0,  k  =  lt  * •  ,  K ,  it  is  a  constrained  least  squares  problem. 


Remark  3.2  Theoretically  the  weighting  masses  can  be  computed  directly  from  the  first  order  reversal 
curves.  This  works  if  the  signals  are  noise-free,  which  is  usually  not  the  case.  Therefore  we  use  the 
least  squares  method. 
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3.1.1  Experimental  results 


In  general  the  magnetostriction  depends  on  both  the  mechanical  pre-stress  a  and  the  magnetic  field  H 
[17,  27].  Pre-stress  is  applied  to  the  magnetostrictive  actuator  through  preloaded  springs  (see  Figure  1) 
and  that  improves  magnetostriction.  The  pre-stress  is  not  adjustable  once  the  actuator  is  manufac¬ 
tured,  and  it  does  not  change  much  during  operation  considering  the  magnitude  of  magnetostriction 
(less  than  1500  parts  per  million  for  Terfenol-D).  Therefore  we  assume  that  the  magnetostriction  is 
only  dependent  on  the  magnetic  field  H,  which  is  determined  by  the  input  current  I. 


For  the  purpose  of  control,  we  define  the  magnetostriction  A  to  be 


A  = 


(5) 


where  lrod  is  the  length  of  the  magnetostrictive  rod  in  the  demagnetized  state,  and  A l  is  the  change 
of  the  rod  length  from  lroci-  Note  that  the  displacement  y  of  the  actuator  head  is  equal  to  A l.  The 
saturation  magnetostriction  As  is  defined  in  an  obvious  way,  which  is  slightly  different  from  that  in 
[28]- 


When  the  input  frequency  is  low,  the  magnetostrictive  hysteresis  is  rate-independent  and  we  can 
relate  y  to  the  bulk  magnetization  M  along  the  rod  direction  by  a  square  law  [4] 


Irod^s  „  t2 

y  =  ~mTm  ’ 


(6) 


and  relate  the  input  current  I  to  the  magnetic  field  H  (assumed  uniform)  along  the  rod  direction  by 


H  =  c0I  +  H{ 


bias  5 


(7) 


where  Ms  is  the  saturation  magnetization,  co  is  the  so  called  coil  factor,  and  Hbias  is  the  bias  field 
produced  by  permanent  magnets  or  a  dc  current.  Hbias  is  necessary  for  generating  bidirectional 
strains.  Hence  we  can  capture  the  hysteretic  relationship  between  y  and  I  by  the  ferromagnetic 
M  —  H  hysteresis.  Venkataraman  employed  a  low  dimensional  ferromagnetic  hysteresis  model  in  [4]. 
Here  we  use  a  Preisach  operator  to  model  M  —  H  hysteresis: 


M  =  T[H,^0\. 


(8) 


Remark  3.3  Due  to  the  thin  rod  geometry,  we  approximate  the  continuum  magnetization  in  the 
magnetostrictive  rod  by  the  bulk  magnetization.  The  square  law  (6)  follows  from  the  continuum  theory 
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DSpace 

ControlDesk 


Figure  7:  Experimental  setup. 

of  micromagnetics,  where  the  magnetoelastic  energy  is  of  the  form  linear  in  the  strain  and  quadratic 
in  the  direction  cosines  of  the  magnetization  vector  [29]. 

Remark  3.4  Mayergoyz  has  shown  that,  the  necessary  and  sufficient  conditions  for  a  hysteretic  non¬ 
linearity  to  be  represented  by  a  Preisach  operator  are  the  wiping-out  property  and  the  congruency 
property  [12].  While  the  wiping-out  property  for  the  ferromagnetic  hysteresis  can  be  directly  verified, 
we  will  indirectly  verify  the  congruency  property  by  a  trajectory  tracking  experiment  based  on  inversion 
of  the  Preisach  operator. 

The  following  parameters  are  available  from  the  manufacturer:  Ms  =  7.87  x  105A/m,  lrod  = 
5.13  x  10~2m,  co  =  1.54  x  104/m.  We  can  easily  identify  As  =  1.313  x  10~3  by  applying  an  input  of 
relatively  large  magnitude.  The  bias  field  Hbias  is  identified  to  be  1.23  x  104A/m. 

Our  experimental  setup  is  as  shown  in  Figure  7  .  DSpace  ControlDesk  is  a  tool  for  real-time 
simulation  and  control.  The  displacement  of  the  actuator  is  measured  with  a  LVDT  sensor,  which  has 
a  precision  of  about  1  ^m. 

The  magnetic  field  H  is  limited  to  [1.57  x  103A/m,  3.25  x  104A/m],  and  we  discretize  the  Preisach 
plane  into  25  levels.  Figure  8  shows  the  distribution  of  the  identified  weighting  masses.  The  constant 
contribution  uq  is  estimated  to  be  4.99  x  105A/m. 

Remark  3.5  Due  to  the  bias  field  H^as  o.nd  the  constraint  on  the  input  current,  we  can  not  trace  the 
major  loop  of  the  M  -  H  hysteresis;  instead  we  can  only  visit  a  certain  region  inside  the  major  loop. 
As  a  result,  the  magnetostrictive  hysteresis  loop  (the  butterfly  curve)  is  asymmetric  (Figure  2). 
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Figure  8:  Distribution  of  the  Preisach  weighting  masses. 

3.2  Inversion  of  the  discretized  Preisach  operator 


General  models  for  smart  actuators  that  capture  both  hysteresis  and  dynamic  behaviors  have  a  cas¬ 
caded  structure,  as  shown  in  Figure  9.  Here  W  is  some  hysteretic  nonlinearity  and  G(s )  represents 
the  transfer  function  of  the  linear  part  in  the  actuator. 


yref 


Controller 


Figure  9:  Controller  design  schematic  [30]. 

A  basic  idea  for  controller  synthesis  for  such  systems  is  to  construct  a  right  inverse  operator  VF-1 
for  W .  Then  v(-)  =  v(-)  and  the  controller  design  problem  is  reduced  to  designing  a  linear  controller 
K(s)  for  the  linear  system  G(s)  (Figure  9). 

The  hysteretic  nonlinearity  W  could  be  a  rate-independent  model,  e.g.,  the  Preisach  operator,  or 
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it  could  be  a  rate-dependent  model,  like  the  dynamic  model  to  be  discussed  in  the  next  section.  In  this 
subsection  and  the  next  subsection,  we  study  inversion  of  the  Preisach  operator.  Preisach  operator  is 
highly  nonlinear,  and  in  general,  we  can  not  find  a  closed-form  formula  for  the  inverse  operator,  unless 
the  density  function  is  of  some  special  form,  as  in  the  work  of  Galinaitis  and  Rogers  [31].  Hughes 
and  Wen  [5,  18]  utilized  the  first  order  reversal  curves  in  computing  the  numerical  inverse  of  the 
Preisach  operator.  This  method  relies  on  measurement  of  all  first  order  reversal  curves  and  involves 
solving  nonlinear  equations.  Natale  and  his  colleagues  proposed  using  another  Preisach  operator  as 
a  “pseudo-compensator”  to  approximate  the  inverse  of  a  Preisach  operator  [10].  The  compensator  is 
“pseudo”  because  it  is  well  known  that  in  general,  the  inverse  of  a  Preisach  operator  is  not  a  Preisach 
operator.  Venkataraman  and  Krishnaprasad  proposed  an  inversion  algorithm  based  on  the  contraction 
mapping  principle  for  a  Lipschitz  continuous  Preisach  operator  [30] . 

In  this  paper  we  exploit  the  piecewise  nronotonicity  property  of  the  Preisach  operator  with  a 
nonnegative  measure,  and  other  structural  properties  inherited  from  the  discretization,  to  develop 
efficient  inversion  algorithms  which  are  implementable  in  real-time.  We  first  study  inversion  of  a 
discretized  Preisach  operator  obtained  as  a  result  of  input  discretization. 

Let  U  be  the  set  of  discrete  input  values,  i.e. ,  U  =  {iq,  1  <  l  <  L  +  1}  with 

.  /  7  i\  r  i  r  'Umax  V"min 

ui  =  umin  +  {l~  1  )ou,  where  du  = - - - . 

jL/ 

Let  Sn  be  the  set  of  input  strings  of  length  n  taking  values  in  U,  i.e.,  if  s  E  Sn,  then  s[i]  6  U,  1  <  i  <  n. 
Let  be  the  set  of  memory  curves  for  the  discretized  Preisach  operator.  Note  that  any  ip  £ 
consists  of  L  segments,  and  each  segment  can  be  either  horizontal  or  vertical  (see  Figure  5(b)  for  an 
example  of  an  element  in  with  L  =  3).  Therefore,  there  are  2L  elements  in 

Inversion  Problem  of  Length  N:  Let  all  measure  masses  be  nonnegative.  Given  an  initial 
memory  curve  ip o  £  'Ld  and  a  desired  output  sequence  z  of  length  N ,  find  s*  £  SN ,  such  that 

max  |r[s*,«l  —  z\i]\  =  min  max  |r[,s, V’olRl  —  z\i]\.  (9) 

1  <i<N  JLJ  seSN  1<1<N  JLJ  LJI 

Remark  3.6  Although  the  Preisach  measure  for  a  discretized  Preisach  operator  is  not  nonsingular, 
the  configuration  of  hysterons  at  any  time  is  fully  captured  by  some  memory  curve  in  and  we  can 
still  write  r [• ,  -0]  for  ip  £ 

Remark  3.7  A  discretized  Preisach  operator  is  not  uonto”  since  its  output  takes  values  in  a  finite 
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set.  Therefore,  we  do  not  seek  an  exact  inverse  in  the  problem  formulation. 


Before  we  present  the  solution  to  the  problem  above,  we  first  look  at  the  case  when  N  =  1: 

Inversion  Problem  of  Length  1:  Let  all  measure  masses  be  nonnegative.  Given  an  initial 
memory  curve  V’o  £  'k d  and  a  desired  output  zq,  find  u*  £  U ,  such  that 

|I>*,^o]  -zo\  =  min|r[u,V’o]  ~  zo\-  (10) 

u€U 

There  is  a  simple  algorithm  for  solving  the  problem  of  length  1,  which  is  based  on  piecewise 
monotonicity  of  the  Preisach  operator  [9] .  We  name  it  the  closest  match  algorithm ,  because  it  always 
generates  an  input  such  that  the  corresponding  output  matches  the  desired  output  most  closely  among 
all  possible  inputs. 

The  idea  of  the  closest  match  algorithm  is  as  follows.  One  can  obtain  the  initial  input  u ^  and  the 
initial  output  from  the  initial  memory  curve  ipo.  Consider  the  case  z^  <  zq  (the  case  z^  >  zq 
is  treated  in  exactly  the  same  way  with  some  obvious  modification).  We  keep  increasing  the  input  by 
one  level  in  each  iteration  until,  say  at  iteration  n,  the  input  u (n')  reaches  umax,  or  the  output  z Fd 
corresponding  to  u ^  exceeds  zq.  For  the  first  case,  the  optimal  input  is  clearly  umax •  For  the  second 
case,  two  candidates  for  the  optimal  input  u*  are  u^n~^  and  u^n\  and  we  take  u*  to  be  the  one  with 
smaller  output  error.  Note  that  we  need  back  up  the  memory  curve  whenever  we  increase  the  input, 
so  that  we  can  always  retrieve  the  consistent  memory  curve  with  u*. 

The  above  algorithm  yields  the  optimal  input  u*  in  at  most  L  iterations.  And  in  each  iteration, 
the  evaluation  of  z ^  is  very  fast  since  the  input  has  changed  by  one  level  and  thus  we  need  only 
update  states  of  hysterons  corresponding  to  that  level.  These  factors  combine  to  make  this  algorithm 
simple  and  efficient. 

The  trajectory  inversion  problem  of  length  N  is  solved  by  combining  the  closest  match  algorithm 
and  the  dynamic  programming  principle. 

Let  5,;  :  'I'j  x  [/  — >  4^  be  the  evolution  map  for  the  memory  curve,  i.e. ,  if  if  £  is  the  initial 
memory  curve,  then  Ed(u,if>)  is  the  new  memory  curve  when  the  input  u  £  U  is  applied. 
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Given  N  and  a  desired  output  sequence  z  of  length  N ,  we  define,  for  1  <  k  <  N, 


k<i<N 


Ml,  s  6  sK-‘+1, 

(11) 

(12) 

where  we  call  Jk  the  cost  function  and  14  the  value  function. 


Proposition  3.1  The  value  functions  14,  1  <  k  <  N,  can  be  solved  successively  via: 

VnW  =  min  \T[u,if]  -  z[N}\,  (13) 

ueu 

Vk (VO  =  min max{ |r[u, ip\  -  z[k\\,Vk+i{Ed(u,i/j))}.  (14) 

ueu 

Define  maps  n(  :  44  —>  U,  1  <  k  <  N,  so  that  vr^(^)  is  the  argrnin  in  (13)  and  (If).  Then  for  the 
inversion  problem  of  length  N ,  ir 1  <  k  <  N,  gives  the  optimal  control  policy  at  time  k. 

Proof.  Straightforward  from  Bellman’s  optimality  principle  [32].  □ 

The  closest  match  algorithm  can  be  used  in  solving  (13)  and  (14).  Proposition  3.1  entails  pre- 
computing  and  storage  of  the  optimal  maps,  which  is  undesirable  when  N  or  the  cardinality  of  44  is 
large.  A  sub-optimal  approach  is  to  decompose  the  inversion  problem  of  length  N  into  N  successive 
inversion  problems  of  length  1,  and  solve  them  using  the  closest  match  algorithm.  Experimental  results 
of  trajectory  tracking  based  on  this  approach  can  be  found  in  [9] . 


3.3  Inversion  of  the  Preisach  operator  with  a  nonsingular  measure 

We  now  discuss  the  inversion  problem  for  a  Preisach  operator  with  a  nonsingular  Preisach  measure. 
In  this  case,  the  Preisach  operator  can  be  inverted  with  arbitrary  accuracy,  and  it  suffices  to  study 
an  inversion  problem  of  length  1:  given  ifo  G  and  M  £  [Mmm,  Mrnax\.  find  H  £  [Hmin,  Hmax\,  such 
that 

M  =  T[H,fi0\, 

where  [Hminj  Hmax\  and  [Mmm ,  Mmax]  are  the  ranges  of  the  input  and  the  output  of  the  Preisach  oper¬ 
ator,  respectively.  The  notation  used  in  this  subsection  is  slightly  different  from  that  in  Subsection  3.2, 
but  it  will  be  consistent  with  the  notation  in  Section  4. 
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Proposition  3.2  Let  the  Preisach  measure  be  nonnegative  and  nonsingular  with  a  density  function 


p.  Let 


A 


u  =  max  {sup 


/  p(/3,  a)d(3,  sup 

po  P 


rot  o 

h 


p(/3,  a)da}  <  oo. 


(15) 


Let  V’o  £  ^  be  the  initial  memory  curve,  and  let  the  input  and  the  output  of  the  Preisach  operator 
corresponding  to  ifo  be  Ho  and  Mq,  respectively.  Given  M  e  \Mmin,  Mmax\,  consider  the  following 
algorithm: 


pj(n+\)  _  jj{n)  _j_ 

M(n+1)  =  Y[H^n+1\lf^} 


(16) 


where  ip^  =  if q,  =  Ho,  M®  =  Mo,  and  if^  is  the  memory  curve  after  {i^(fc)}”=1  is  applied. 
Then  M ^  — >  M  as  n  — ►  oo. 


Proof.  The  proposition  follows  directly  from  piecewise  monotonicity  and  continuity  of  the  Preisach 
operator.  □ 


Remark  3.8  The  algorithm  (16)  also  appeared  in  [30],  where  approximate  inversion  of  the  Preisach 
operator  was  studied  for  the  class  of  continuous,  piecewise  monotone  functions. 

What  we  have  identified  in  Subsection  3.1  is  a  set  of  Preisach  weighting  masses  and  the  corre¬ 
sponding  Preisach  measure  is  not  nonsingular.  We  can  obtain  a  nonsingular  Preisach  measure  vp  by 
assuming  that  each  identified  mass  is  distributed  uniformly  over  the  corresponding  cell  in  the  dis¬ 
cretization  grid.  Note  that  the  diagonal  cells  are  triangular,  while  other  cells  are  square  (refer  to 
Figure  5(a)).  The  density  function  pp  corresponding  to  vp  is  piecewise  uniform,  which  enables  us  to 
solve  the  inversion  problem  exactly  (in  a  finite  number  of  iterations),  as  described  next. 

We  consider  the  case  M  >  Mo  and  the  other  case  can  be  treated  analogously.  It’s  obvious  that 

-  /n\ 

H  >  Hq  and  we  will  increase  the  input  in  every  iteration.  At  iteration  n,  let  d\  >0  be  such  that 

H^l  +  d![L'>  equals  the  next  discrete  input  level,  and  let  >  0  be  the  minimum  amount  such  that 

applying  +  dP1'1  will  eliminate  the  next  corner  of  the  memory  curve  (see  Figure  10  for  illustration). 

Since  pp  is  piecewise  constant,  for  d  <  min{d(jn\  dij1'1},  we  have 

T[H^  +  d,if{n)]  -  T[H^n\  ^(n)]  =  4n)d2  +  4n)d, 
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Figure  10:  Illustration  of  and  d1^1  (L  =  8). 

(n)  ( n ) 

where  a)  ,  a^'  —  0  can  be  computed  from  fip,  and  the  square  term  is  due  to  the  contribution  from 

(n) 

the  triangular  region  inside  the  diagonal  cell.  Let  dy0  '  be  such  that 


M) 


m  -  y[h^\^)  =  4n)(<4n))2  +  <4n)4n). 


The  inversion  algorithm  now  works  as  follows: 

d(n)  =  min{dg!\  d^'\  d^} 

<  H^n+V>  =  H^>  +  dW  •  (17) 

M(n+1)  =  r[H^n+1\'lp^} 

If  at  iteration  n*,  d^1*^  =  d^1  \  then  the  iteration  stops  and  H  =  iL(n*+1).  Let  nc('i/,o)  be  the  number 
of  corners  of  i(jq,  and  L  be  the  discretization  level  of  the  Preisach  plane.  Then  it  is  easy  to  see  the 
algorithm  (17)  yields  the  (exact)  solution  in  no  more  than  n  =  nc(tpo)  +  L  iterations. 

Figure  11  shows  the  result  of  an  open-loop  tracking  experiment  using  the  algorithm  (17).  The 
desired  trajectory  was  obtained  from  the  output  of  a  Van  der  Pol  oscillator  to  make  the  tracking 
task  challenging.  In  Figure  11,  the  displacement  trajectories  (both  the  desired  and  the  measured), 
the  tracking  error,  and  the  input  current  computed  based  on  the  inversion  algorithm  are  displayed. 
The  overall  performance  is  satisfactory  since  the  error  magnitude  is  less  than  3  /um  most  of  the  time 
with  a  tracking  range  of  60  /im.  We  can  see  that  the  tracking  error  slightly  exceeds  4  /x nr  when  the 
desired  output  (and  thus  the  input)  undergoes  abrupt  changes,  in  which  case  the  rate-independence 
assumption  no  longer  holds. 

The  inversion  problem,  as  discussed  here  and  in  the  literature,  is  a  trajectory  inversion  problem, 
since  we  want  to  find  an  input  sequence  such  that  the  output  sequence  matches  (most  closely)  the 
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Figure  11:  Trajectory  tracking  based  on  inversion  of  the  Preisach  operator. 

desired  one.  In  some  applications,  e.g.,  micro-positioning,  one  is  more  interested  in  finding  an  input 
sequence  such  that  the  final  value  of  the  output  sequence  matches  (most  closesly)  the  desired  value. 
This  problem,  called  the  value  inversion  problem,  was  formulated  and  studied  in  [24]. 


4  Identification  and  Inversion  of  a  Dynamic  Hysteresis  Model 

4.1  A  dynamic  magnetostrictive  hysteresis  model 

When  the  input  frequency  gets  high,  the  magnetostrictive  hysteresis  is  rate-dependent  (Figure  12)  due 
to  the  eddy  current  losses  and  the  magnetoelastic  dynamics  of  the  magnetostrictive  rod.  Venkataraman 
and  Krishnaprasad  proposed  a  bulk  hysteresis  model  for  a  thin  magnetostrictive  actuator  based  on 
energy  balancing  principles  [33,  4].  The  model  has  a  cascaded  structure  as  shown  in  Figure  13.  W 
takes  care  of  the  M  -  H  hysteresis  and  the  eddy  current  losses,  and  G(s )  is  a  second  order  linear 
system  modelling  the  magnetoelastic  dynamics. 
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Figure  12:  The  rate-dependent  magnetostrictive  hysteresis. 


Figure  13:  Model  structure  of  a  magnetostrictive  actuator. 


We  now  have  a  closer  look  at  the  block  W.  Due  to  the  finite  resistivity  of  the  magnetostrictive 
material,  there  are  eddy  currents  circulating  inside  the  rod.  One  way  to  represent  the  eddy  current 
losses  is  to  place  a  resistor  Reddy  hr  parallel  with  a  hysteretic  inductor  [28,  4],  as  shown  in  Figure  14. 
We  note  that  this  is  a  phenomenological  approach  and  the  underlying  details  of  the  eddy  current 
dynamics  are  ignored  here.  We  assume  that  the  magnetic  flux  density  B  is  uniform  over  the  cross 
section  of  the  magnetostrictive  rod.  Then  the  voltage  V  across  the  nonlinear  inductor  is  NmAm^§-, 
where  Nm  is  the  number  of  turns  of  the  coil,  and  Am  is  the  cross  sectional  area  of  the  rod.  Let  /  be 
the  input  current  applied  to  the  actuator,  and  I\  be  the  current  flowing  in  the  inductor  branch.  Since 


V  =  (I  -  h) Reddy,  we  have 


dB_ 

dt 


Reddy  ,T  T  \ 
N  A  1 

1  Vm^1m 


(18) 


In  SI  units,  B  =  +  M),  where  /io  =  47r  x  10  'Henry/m  is  the  permeability  of  vacuum.  H  is 

related  to  I\  via  H  =  cq/i  +  H^ias. 
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Figure  14:  Representation  of  eddy  current  losses  in  a  magnetostrictive  actuator  [4]. 

Remark  4.1  Again,  considering  the  thin  rod  geometry  of  the  actuator,  the  continuum  magnetization 
and  the  magnetic  flux  density  are  approximated  by  their  bulk  values  along  the  rod  direction,  and  the 
magnetoelastic  dynamics  is  lumped  into  a  second  order  linear  system  G(s) .  Since  the  goal  is  to  develop 
a  phenomenological  model  suitable  for  real-time  control  applications,  the  approximation  is  necessary. 

The  constitutive  relationship  between  M  and  H  was  modeled  by  a  low  dimensional  bulk  ferro¬ 
magnetic  hysteresis  model  in  [4]  and  that  led  to  an  overall  model  described  by  switching  ordinary 
differential  equations.  We  use  a  Preisach  operator  T  with  a  nonnegative  Preisach  measure  to  model 
M  —  H  hysteresis  and  obtain  the  following  new  dynamic  model  for  the  block  W : 


H(t)  +  M(t)  =  ci  (1(f) 
M(t)  =T[H(-)^0](t) 


where  ipo  represents  the  initial  memory  curve  and 


H(t)  —  Hbias  \ 

co  1 


Cl 


A 


R, 


eddy 


VoNmAr 


(19) 


G(s)  has  a  state  space  representation  [33,  4]  (after  some  manipulations): 

y(t)  +  2 &oy(t)  +  y(t)  =  u°l^Xs  M2(t),  (20) 

where  y  is  the  displacement,  u>o  =  27t/o,  /o  is  the  first  resonant  frequency  of  the  actuator,  £  is  the 
damping  coefficient,  lrod  is  the  length  of  the  rod,  As  is  the  saturation  magnetostriction  and  Ms  is  the 
saturation  magnetization. 
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Figure  15:  The  block  diagram  of  the  block  W. 

We  note  that  the  model  (19)  and  (20)  degenerates  to  the  rate-independent  model  in  Section  3  if 
we  set  the  derivatives  in  these  equations  to  0. 

In  (19),  the  Preisach  operator  is  coupled  to  an  ordinary  differential  equation  (ODE)  in  an  unusual 
way.  Eq.  (19)  can  be  viewed  as  a  special  nonlinear  feedback  system  (Figure  15),  and  it  presents 
interesting  problems  in  analysis  and  computation  [24]. 

The  following  theorem  shows  that  the  model  (19)  is  well-posed.  The  proof  can  be  found  in  [24,  22]. 

Theorem  4.1  If  the  Preisach  measure  v  is  nonnegative  and  nonsingular,  and  /(•)  is  piecewise  contin¬ 
uous,  then  for  any  initial  memory  curve  tpo,  for  any  T  >  0,  there  exists  a  unique  pair  {H(-),M(-)}  £ 
C([0,T])  x  C([0,T])  satisfying  (19)  almost  everywhere. 

4.2  Identification  of  the  dynamic  model 

We  can  identify  a  discrete  approximation  to  the  Preisach  measure  using  the  scheme  presented  in 
Subsection  3.1.  Furthermore,  we  can  obtain  a  nonsingular  Preisach  measure  from  the  collection  of 
identified  weighting  masses  by  distributing  each  mass  uniformly  over  the  corresponding  cell  in  the 
discretization  grid,  as  discussed  in  Subsection  3.3.  This  is  important  since  the  well-posedness  of  (19) 
requires  the  Preisach  measure  to  be  nonsingular. 

The  following  parameters,  in  addition  to  those  listed  in  Subsection  3.1,  are  available  from  the 
manufacturer:  Nm  =  1300,  Am  =  2.83  x  10-5m2.  To  estimate  the  first  resonant  frequency,  we 
apply  sinusoidal  inputs  of  the  same  amplitude  but  different  frequencies  and  measure  the  amplitudes 
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Figure  16:  Displacement  amplitude  vs.  input  frequency. 

of  the  displacement.  Figure  16  displays  the  displacement  amplitudes  at  different  frequencies  and  we 
determine  the  first  resonant  frequency  to  be  392  Hz. 

The  most  difficult  parameters  to  identify  are  Reddy  and  £  due  to  the  coupling  of  (19)  and  (20).  A 
theoretical  value  of  Reddy  can  be  computed  with  the  formula  Reddy  =  ^  [4],  where  p  is  the 

resistivity  of  the  magnetostrictive  material,  b  and  a  are  the  outer  and  inner  radii  of  the  magnetostrictive 
rod.  We  obtain  an  upper  bound  R  of  Reddy  by  letting  a  =  0.  We  then  discretize  [0,  R]  and  denote  the 
mesh  points  by  R\ddy>  *  =  "  >  Ah  The  discretization  need  not  be  uniform  and  we  make  it  finer  in 

the  region  where  the  dynamics  of  (19)  is  more  sensitive  to  Reddy- 

We  observe  a  periodic  motion  of  the  actuator  head  when  a  periodic  input  is  applied.  The  existence 
of  periodic  solutions  of  the  model  under  periodic  forcing  has  been  proved  in  [24] .  In  addition,  numerical 
simulation  shows  that  the  steady-state  solution  of  (19)  and  (20)  is  periodic  when  /(•)  is.  These 
observations  motivate  the  following  scheme  to  identify  Reddy  and  £: 

•  Step  1.  We  apply  a  sinusoidal  current  (with  some  dc  shift)  /(•)  with  frequency  /  to  the  actuator 
and  measure  the  phase  lag  6yj  between  the  fundamental  frequency  component  of  y(-)  and  /(•); 

•  Step  2.  For  each  R^ddy,  we  numerically  integrate  (19)  with  /(•)  as  the  input,  and  calculate  the 
phase  lag  9^2  j  between  the  fundamental  frequency  component  of  M2(-)  and  /(•). 

•  Step  3.  The  difference  9yj  —  0M2d  is  considered  to  be  the  phase  lag  between  the  fundamental 
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frequency  component  of  y(-)  and  that  of  M2(-)  in  (20),  from  which  we  can  compute  £®. 


Remark  4.2  The  idea  of  relating  the  phase  shift  between  the  output  and  the  input  to  hysteresis  can 
also  be  found  in  [34]-  We  note  that  in  general,  the  phase  lag  depends  highly  nonlinearly  on  the  initial 
conditions  (and  the  amplitude  and  frequency  of  I{-)),  so  we  should  make  sure  that  the  initial  conditions 
in  numerical  simulation  are  consistent  with  experiment  conditions. 


We  repeat  the  above  procedure  (Step  1  to  Step  3)  K  times  with  different  input  frequencies  and 
denote  the  damping  coefficients  as  {^*'*}(L1  for  R^}dy-  If  R^ddy  is  close  to  the  true  parameter  Reddy, 
&  should  not  vary  much  with  k.  Therefore,  we  pick  i*  e  {l,---  ,N}  such  that  has  the 

minimum  variance,  and  estimate  Reddy  via  Reddy  =  Reddy  and  ^  £  be  the  mean  of  {£*.*  ^}k=v 

Figure  17  shows  the  variation  of  £  with  respect  to  frequency  for  different  RKeddy  s.  The  parameters 
are  determined  to  be  Reddy  =  7017,  £  =  0.7783.  Figure  18  compares  the  rate-dependent  hysteresis  loops 
measured  in  experiments  to  those  obtained  through  simulation  based  on  the  identified  parameters. 
We  see  that  the  simulation  results  agree  with  the  experimental  results  reasonably  well  up  to  200  Hz. 
Beyond  200  Hz,  they  still  qualitatively  agree  although  they  match  worse.  This  reveals  that  further 
details  of  the  eddy  currents  and  the  magnetoelastic  dynamics  need  to  be  brought  into  the  picture  to 
fully  capture  the  dynamic  behaviors  at  very  high  frequency. 


Figure  17:  Identification  of  Reddy  and  £. 
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Figure  18:  Model  validation.  Solid  line:  experimental  measurement;  Dashed  line:  numerical  prediction 
based  on  the  identified  parameters. 

4.3  Inversion  of  the  dynamic  model 


Given  an  initial  memory  curve  i/jo  G  'h  and  a  desired  trajectory  M(-)  (obtained,  e.g.,  as  an  output  of 
K(s)  in  Figure  9),  the  inversion  problem  for  (19)  is  to  find  /(•),  such  that  the  corresponding  output 
of  W  is  M(-). 


We  propose  the  following  (formal)  inverse  scheme  for  (19): 

H(t)  =  r~1[M(.)^0}(t) 


m  =  ±(H(t)+M(t))  + 


co 


(21) 


Due  to  the  uniqueness  of  solution  to  (19),  we  expect  the  output  M(-)  under  /(•)  to  agree  with  M(-). 

In  implementation,  we  can  use  the  schemes  discussed  in  Section  3  for  inversion  of  the  Preisach 
operator  (the  first  equation  in  (21)),  and  M  and  H  in  the  second  equation  in  (21)  are  approximated 
by  the  finite  difference  method. 
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Trajectory  tracking  experiments  have  been  carried  out  to  validate  the  model  and  examine  the 
performance  of  the  inverse  scheme  (21).  To  highlight  the  idea  of  inverse  control,  we  have  picked 
£/(•)  from  the  space  of  attainable  output  trajectories  y{  )  of  G(s)  (recall  (20))  under  some  control 
'«(•)  €  C([0,r])  with  0  <  u(t)  <  Mg,  Vi  G  [0, T],  In  this  case,  M  can  be  directly  computed  as:  Vi, 
M(t)  =  y/u(t),  where 

M2 

u{t)  =  2  s  (§(t)  +  2&0§(t)  +  uly(t)). 

We  compare  the  inverse  control  scheme  (21)  based  on  the  dynamic  model  with  a  proportional 
feedback  scheme  and  two  other  inverse  schemes.  For  the  other  two  inverse  schemes,  the  first  one  is 
based  on  the  static  hysteresis  model  in  Section  3,  and  the  second  one  is  based  on  a  non-hysteretic 
model  where  the  input-output  relationship  is  approximated  by  a  single-valued  function 

y(t)  =  — 7.44/3(i)  -  2.63/2(i)  +  40.811(f)  +  30.34.  (22) 

The  non-hysteretic  model  (22)  was  obtained  through  the  least  squares  method. 

Experimental  results  are  shown  in  Figures  19-22.  The  sampling  period  used  was  5  x  10-5  second. 
In  each  figure,  the  displacement  trajectories  (both  the  desired  and  the  measured),  the  tracking  error, 
and  the  input  current  applied  are  displayed.  We  can  see  that  the  performance  of  the  inverse  control 
scheme  based  on  the  dynamic  hysteresis  model  is  very  satisfactory  (Figure  19).  Although  the  controller 
parameter  has  been  carefully  tuned,  the  performance  of  the  pure  feedback  scheme  is  poor  (Figure  20). 
This  highlights  the  need  for  hysteresis  compensation.  The  other  two  inverse  schemes  also  perform 
worse  than  the  scheme  (21),  as  one  can  see  in  Figure  21  and  Figure  22. 
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Figure  19:  Experimental  results  of  trajectory  tracking  using  inverse  control  based  on  the  dynamic 
hysteresis  model. 


Figure  20:  Experimental  results  of  trajectory  tracking  using  proportional  feedback  control. 
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Figure  21:  Experimental  results  of  trajectory  tracking  using  inverse  control  based  on  the  static  hys¬ 
teresis  model. 


Figure  22:  Experimental  results  of  trajectory  tracking  using  inverse  control  based  on  the  non-hysteretic 
model. 
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5  Conclusions 


In  this  paper  we  studied  modeling,  identification  and  control  of  hysteresis  in  smart  actuators  through 
the  example  of  controlling  a  commercially  available  magnetostrictive  actuator. 

Hysteresis  in  a  number  of  smart  materials  can  be  modeled  by  the  Preisach  operator.  We  devel¬ 
oped  fast  inversion  algorithms  for  the  Preisach  operator  by  exploiting  structures  and  properties  of  the 
Preisach  operator.  We  also  presented  a  novel  dynamic  model  for  the  rate-dependent  magnetostrictive 
hysteresis.  This  model  provides  insight  into  modeling  of  rate-dependent  hysteresis  in  other  smart  ma¬ 
terials  in  the  frequency  ranges  of  practical  interest.  Parameter  identification  and  inverse  compensation 
for  the  dynamic  model  were  also  discussed  in  detail.  Extensive  experimental  results  have  shown  that 
the  model  we  proposed  can  capture  the  hysteretic  and  dynamic  behaviors  of  the  actuator,  and  that 
our  identification  and  inverse  control  schemes  are  effective.  It  is  noteworthy  that  we  achieved  almost 
perfect  tracking  of  an  irregular  and  relatively  high  frequency  signal  for  the  full  operating  range  of  the 
actuator. 

The  approaches  presented  in  this  paper  are  very  general  and  they  are  applicable  to  control  of  a 
wide  class  of  smart  actuators. 

Due  to  the  open  loop  nature  of  the  inverse  control  scheme,  its  performance  is  susceptible  to 
model  uncertainties  and  to  errors  introduced  by  the  inversion  schemes.  To  address  this  problem,  we 
have  combined  inverse  compensation  with  the  li  robust  control  theory  and  proposed  a  robust  control 
framework  for  smart  actuators.  This  will  be  reported  in  [2]. 
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